/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name : f_pow.c
 *  Function  :
 *        f_modef function returns powed value whose original value is given.
 *        This file is copied from math.lib of 
 *                  CC33 tool(CC33v40.exe + CCPA45.exe, 
 *                  math of ansilib33v25 )
 *
 *  original file's Revision  :
 *     1993/08/31    first release                 K.Shinozuka
 *     1993/12/17    change include file           K.Shinozuka
 *     1997/03/17     modified                     T.Katahira
 *                      change function declaration
 *                                    to the ANSI style
 *     1997/03/25    change static to const static      M.Kudo
 *  Revision  :
 *      2003/ 3/12    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/
#ifdef UNIX
#include "f_math.h"
#include "f_smcvals.h"

#else
#include <f_math.h>
#include <f_smcvals.h>
#endif

#define CHECKY 100.0f

static float f_Int( float, int * );


float
f_pow( float x, float y )

/********************************************************************/
/*                                                                  */
/*  The pow function calculates x raised to the power y.            */
/*                                                                  */
/*  RETURN  normal: x raised to the power y                         */
/*          error : NaN (input error)                               */
/*                  - x is negative and y is not an integral valuie */
/*                  - x is zero and y is less than or equal to zero */
/*                  Inf, -Inf ( overflow )                          */
/*                                                                  */
/*     This function uses the external functions,                   */
/*            exp.c, log.c , ldexp.c.                               */
/*                                                                  */
/********************************************************************/

/* ALGORITHM
 * compute pow(x,y) = x^y as exp( y * log(x) )
 * 1.check the sign of pow(x,y)
 *     x is negative and y is odd number --- pow(x,y) is negative
 *                                           let x = -x and sign = -1
 *     otherwise, pow(x,y) is positive and sign = 1
 * 2.|y| >= CHECKY
 *     return( sign * exp( y * log( x ) )
 * 3.let x to be 2^xexp * xfrac  ( sqrt(2)/2 <= xfrac < sqrt(2) )
 *     pow(x,y) = exp( y * xexp * ln2 ) * exp( y * log( xfrac ) )
 * 4.separate y*xexp into integral value and fraction
 *     y * xexp = i_yx + f_yx    ( 0 <= f_yx < 1 )
 * 5.breake y into integral value and fraction
 *     y = i_y + f_y    ( 0 <= f_y < 1 )
 * 6.pow(x,y) = sign * exp( i_yx * ln2 + i_y * log( xfrac ) ) 
 *                        *  exp( f_yx * ln2 + f_y * log( xfrac ) )
 *            = sign * 2^i_yx * xfrac^i_y * exp( f_yx * ln2 * f_y * log( xfrac ) ) 
 */

/*  the architecture of float floating point
 *
 *   0 30         23 22                0  32 bit
 *   -----------------------------------
 *  | |   exponent  |  fraction         |
 *   -----------------------------------
 *
 *  |               |                   |
 *  |    8 bits     |    23 bits        |
 *
 *         bit    31         sign bit        (  1 bit  )
 *             30 - 23      exponent part    (  8 bits )
 *             22 -  0      fraction part    ( 23 bits )
 *
 */


	/* double          x;           the base                                     */
	/* double          y;           the power                                    */
	{
	 FLT_LNG        ckarg;       /* exchange the type to check the argument      */ 
	 FLT_LNG        power;       /* pow( x, y )                                  */
	 float          i_y;         /* intetger part of y                           */
	 float          f_y;         /* fraction part of y                           */
	 float          xfrac;       /* sqrt(2)/2 <= xfrax < sqrt(2)                 */
	 int            xexp;        /* x = 2^xexp * xfrac                           */
	 float          p_xfrac;     /* xfrac^i_y                                    */
	 float          i_yx;        /* integer part of y * xexp                     */
	 float          exp2;        /* exp( frac(y*xexp) * ln2 + f_y * log( xfrac ) */
	 register int   sign;        /* sign of x, 1 ( x>0 ) or -1 ( x<0 )           */
	 int            odd;         /* y is odd number odd = -1; else odd = 1       */
	 register short shorty;      /* change float y to short int y                */
	 const static float  ln2 = 0.69314718055994530942f;

/* check the argument */	
	 ckarg._F = x;
	 F_CHECK_ARG( &ckarg );                        /* check NaN */
	                                               /* error: error code 33 domain error */   
	 if ( ckarg.st._LL == f_P_INF.st._LL ) {       /* x = Inf */
	    ckarg._F = y;
	    F_CHECK_ARG( &ckarg );                     /* check NaN */
	                                               /* error: error code 33 domain error */
	    if ( ckarg.st._LL == f_P_INF.st._LL ) {    /* y = Inf */
	       errno = ERANGE;                         /* 34: range error */
	       return( f_P_INF._F );
	    } else if ( ckarg.st._LL == f_N_INF.st._LL ) {  /* y = -Inf */
	       return( 0.0f );
	    } else if ( y > 0.0f ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( f_P_INF._F );
	    } else if ( y < 0.0f ) {
	       return( 0.0f );
	    } else {
	       return( 1.0f );                         /* INF^0   */
	    }
	 } else if ( ckarg.st._LL == f_N_INF.st._LL ) { /* x = -Inf */
	    ckarg._F = y;
	    F_CHECK_ARG( &ckarg );                     /* check NaN */
	                                               /* error: error code 33 domain error */
	    if ( ckarg.st._LL == f_P_INF.st._LL ) {    /* y = Inf */
	       errno = EDOM;                           /* 33: domain error */
	       return( f_NAN._F );
	    } else if ( ckarg.st._LL == f_N_INF.st._LL ) {  /* y = -Inf */
	       return( 0.0f );
	    } else if ( (i_y = f_Int( y, &odd )) != y ) {  /* y is not integer */
	       errno = EDOM;                           /* 33: domain error */
	       return( f_NAN._F );
	    } else if ( y > 0.0f ) {
	       if ( odd == 1 ) {                       /* y is even number */
	          errno = ERANGE;                      /* 34: range error */
	          return( f_P_INF._F );
	       } else {
	          errno = ERANGE;                      /* 34: range error */
	          return( f_N_INF._F );
	       }
	    } else if ( y < 0.0f ) {
	       return( 0.0f );
	    } else {
	       return( 1.0f );
	    }
	 }

	 ckarg._F = y;
	 F_CHECK_ARG( &ckarg );                        /* check NaN */
	 if ( ckarg.st._LL == f_P_INF.st._LL ) {       /* y = Inf */
	    if ( x <= -1.0f ) {
	       errno = EDOM;                           /* 33: domain error */
	       return( f_NAN._F );
	    } else if ( x > 1.0f ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( f_P_INF._F );
	    } else if ( x != 1.0f ) {
	       return( 0.0f );
	    } else {
	       return( 1.0f );
	    }
	 } else if ( ckarg.st._LL == f_N_INF.st._LL ){ /* y = -Inf */
	    if ( ( x < -1.0f ) || ( x > 1.0f ) ) {
	       return( 0.0f );
	    } else if ( ( x >= -1.0f ) && ( x < 0.0f ) ) {
	       errno = EDOM;                           /* 33: domain error */
	       return( f_NAN._F );
	    } else if ( ( x >= 0.0f ) && ( x < 1.0f ) ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( f_P_INF._F );
	    } else {
	       return( 1.0f );
	    }
	 }

	 if ( x == 0.0f ) {                       /* zero */
	    if ( y <= 0.0f ) {                    /* y is less than or equal to zero */
	       errno = EDOM;                      /* 33: domain error */
	       return( f_NAN._F );
	    } else {                              /* y is positive */
	       return( 0.0f );
	    }
	 }

	 if ( x == 1.0f ) {
	    return( 1.0f );
	 } else if( x > 0.0f ) {                  /* x is positive */
	    sign = 1;
	 } else if ( (i_y = f_Int( y, &odd )) != y ) {   /* y is not integer */
	    errno = EDOM;                         /* 33: domain error */
	                                          /* avoid imaginary part */
	    return( f_NAN._F );
	 } else {                                 /* y is integer */
	    x = -x;
	    sign = odd;
	 }
	 if ( y == 0.0f ) {
	    return( 1.0f );
	 }

	 if ( y > CHECKY || y < - CHECKY ) {
	    power._F  = f_exp( y * f_log( x ) );
	    if ( power.st._LL > f_P_INF.st._LL ) {  /* overflow in exp function */
	       errno = ERANGE;                    /* 34: overflow */
	       if ( sign == -1 ) {                /* negative overflow */
	          return( f_N_INF._F );
	       } else {
	          return( f_P_INF._F );
	       }
	    } else {
	    
	       if ( sign == -1 ) 
	          return( -power._F );
	       else
	          return( power._F );
	    }
	 }

	 exp2 = 0.0;
	 i_yx = 0.0f;
	xfrac = f_frexp( x, &xexp );
	 if ( xexp != 0 ) {
	    i_yx = f_Int( y * (float)xexp, &odd );   /* get integer */
	    exp2 = ( y * (float)xexp - i_yx ) * ln2; // get faraction //
	 }

/* calculate sin * 2^i_yx * p_xfrac * exp( exp2 ) */
	 p_xfrac = 1.0f;
	 if ( xfrac != 1.0f ) {                    /* calculate 2^shorty, exp2 */
	    shorty = (short) y;
	    f_y = y - (float)shorty;
	    if ( f_y != 0.0f ) {
	       exp2 = exp2 + f_y * f_log( xfrac );
	    }

	    if ( shorty < 0 ) {                   /* y is negative number */
	       shorty = -shorty;
	    }
	    while ( shorty != 0 ) {               /* y is not zero */
	       if( shorty & 1 ) {                 /* the least significant bit is HIGH */
	          p_xfrac = p_xfrac * xfrac;      /* get xfrac^shorty */
	       }
	       xfrac = xfrac * xfrac;
	       shorty = ( shorty >> 1 );  
	    }
	    if ( y < 0.0f ) {
	       p_xfrac = 1 / p_xfrac;             /* xfrac^(-shorty) = 1/xfrac^shorty */
	    }
	 }

	 p_xfrac = p_xfrac * f_exp( exp2 );
	 power._F = f_ldexp( p_xfrac, (int)i_yx );  /* 2^i_yx * p_xfrac */
	 if ( power.st._LL > f_P_INF.st._LL ) {   /* overflow in exp function */
	    errno = ERANGE;                       /* 34: overflow */
	    if ( sign == -1 ) {                  /* negative overflow */
	       return( f_P_INF._F );
	    } else {
	       return( f_N_INF._F );
	    }
	 }

	 // call other instead of cast, mutiply emulib
	 //return( sign * power._F );
	 if ( sign == -1 ) {
	   return( -power._F );
	 } else {
	   return( power._F );
	 }


}


static float
f_Int( x, odd )

/************************************************************************/
/*                                                                      */
/*  The _Int function computes the largest integral value not greater   */
/*  than x.  If x is less than 1, odd code is 0. If return value is     */
/*  odd, the code is -1, otherwise the code is 1. The code is stored in */
/*  the object pointed by odd.                                          */
/*  The argument x is positive number, not zero, Inf.                   */
/*                                                                      */
/*  RETURN    the largest integral value not greater than  x            */ 
/*                                                                      */
/************************************************************************/

/* ALGORITHM
 * 1.separate x to fraction part and exponent part
 *     x = 2^xexp * xfrac
 *       xfrac: [1,2)
 * 2.if xexp is more then 51, x doesn't have fractional bits
 * 3.if xexp is less than 0, x doesn't have integral bits
 * 4.else, x has both fractional and integral bits
 * 5.clear fractional bits
 * 6.check the lowest bit of integral part, and set the odd code
 */

	float          x;           /* positive and not Inf number         */
	int            *odd;        /* the pointer, the odd code is stored( when even, *odd = 1, else *odd = -1) */
	{
	 FLT_LNG       xfrac;       /* fraction part of x                  */
	 FLT_LNG       ipart;       /* integral number                     */
	 register long xexp;        /* exponent part of x                  */

	 xfrac._F = x;                
	 xexp = (( xfrac.st._LL & F_MASK_EXP ) >> F_L_FRAC );    /* get exponent part  */
	                                          /* MASK_EXP = 0x7ff00000 */
	 xexp = ( xexp - ((unsigned long)F_NUM1_L >> F_L_FRAC) );    /* get off offset */
	                                          /* NUM1_L = 0x3ff00000 */
	 xfrac.st._LL = ( xfrac.st._LL &  F_MASK_FRAC );    /* get fraction part  */
	                                          /* MASK_FRAC = 0x000fffff */
	 if ( xexp >= F_N_FRAC ) {                  /* no fractional bits  */
	    ipart._F = x;                         /* set integer */
	    if ( xexp == F_N_FRAC ) {               /* the least significant bit is 2^0 */
	       if ( ( xfrac.st._LL & 1 ) == 0 ) { /* the least significant bit is LOW */
	          *odd = 1;
	       } else {
	          *odd = -1;
	       }
	    } else {
	       *odd = 1;
	    }
	 } else if ( xexp < 0 ) {                 /* x < 1 : all fractional bits */
	    ipart._F = 0.0f;                       /* set integer */
	    *odd = 0;

	 } else {
	    ipart._F = x;

	    ipart.st._LL = ( ipart.st._LL & ~( (unsigned long)F_MASK_FRAC >> xexp ));  /* clrear fractional bits */ 

	    if ( (( ipart.st._LL & ((unsigned long)F_MASK_FRAC >> (xexp -1))) == 0 ) && ( xexp != 0 ) ) { 
	                                          /* 2^0 is LOW */
	       *odd = 1;
	    } else if ( xexp != 0 ) {
	       *odd = -1;
	    } else {                           /* return value is 1.0 */
	       *odd = -1;
	    }
	 }

	 return( ipart._F );
	}
